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We analytically derive mean-field models for all-to-all coupled networks of heterogeneous, 
adapting, two-dimensional integrate and fire neurons. The class of models we consider 
includes the Izhikevich, adaptive exponential and quartic integrate and fire models. The 
heterogeneity in the parameters leads to different moment closure assumptions that can 
be made in the derivation of the mean-field model from the population density equation 
for the large network. Three different moment closure assumptions lead to three different 
mean-field systems. These systems can be used for distinct purposes such as bifurcation 
analysis of the large networks, prediction of steady state firing rate distributions, 
parameter estimation for actual neurons and faster exploration of the parameter space. We 
use the mean-field systems to analyze adaptation induced bursting under realistic sources 
of heterogeneity in multiple parameters. Our analysis demonstrates that the presence of 
heterogeneity causes the Hopf bifurcation associated with the emergence of bursting to 
change from sub-critical to super-critical. This is confirmed with numerical simulations of 
the full network for biologically reasonable parameter values. This change decreases the 
plausibility of adaptation being the cause of bursting in hippocampal area CA3, an area 
with a sizable population of heavily coupled, strongly adapting neurons. 

Keywords: integrate-and-fire neuron, mean-field model, hippocampus, bifurcation analysis, bursting 



1. INTRODUCTION 

As computers become more powerful, there is a move to numer- 
ically simulate larger and larger model networks of neurons 
(Izhikevich and Edelman, 2008). While simulation is useful for 
confirming observed behavior it is not as helpful in determining 
the mechanisms underlying the behavior. The tools of dynami- 
cal systems theory, such as bifurcation analysis can be useful in 
this regard when studying single neuron or small network mod- 
els. However, they are not viable for large networks, especially if 
the neurons are not identical. Thus, a common approach is to 
try to extrapolate large network behavior from detailed analysis 
of the behavior of individual cells or small networks (Skinner et 
al., 2005). This can be problematic as networks can have behav- 
ior that is not present in individual cells. For example, individual 
neurons that are only capable of tonic firing when isolated may 
burst when coupled in a network (van Vreeswijk and Hansel, 
2001). Further, large networks may exhibit behavior not present 
in smaller networks. For example, Dur-e-Ahmad et al. studied 
bursting in networks ranging in size from 2 cells to 100. They 
found that bursting occurred in a larger range of parameters for 
larger networks (Dur-e-Ahmad et al., 2012, Figure 7). 

Given the role of bursts in networks of neurons, it is impor- 
tant to understand how a network transitions (bifurcates) from a 
non-bursting behavior, to bursting. Bursting has been suggested 
to be a fairly important and information dense firing mode for 
neurons. For example, single bursts can induce long term poten- 
tiation and depression in the hippocampus, which are important 
processes for learning and memory (Lisman, 1997). Additionally, 
bursts have been found to carry more information about an 



animal's position in space than isolated spikes alone as place fields 
have been found to be more accurately defined when consider- 
ing bursts alone (Lisman, 1997). Additionally, the mechanism we 
analyze in this paper, adaptation induced bursting has also been 
suggested as a biologically plausible mechanism for the generation 
of grid cells via oscillatory interference of the bursting (Zilli and 
Hasselmo, 2010). However, adaptation induced bursting is a net- 
work level phenomenon and we cannot apply bifurcation analysis 
directly to a network of adapting neurons. 

Mean-field theory offers one approach to bridge this gap. 
In applying this theory, one usually derives (or suggests) a low 
dimensional system of differential equations that govern the 
moments of different variables across the network (Bressloff, 
2012). For example, for a network of all-to-all coupled Izhikevich 
neurons (Izhikevich, 2003), one can derive a two dimensional sys- 
tem of differential equations for the mean adaptation variable and 
the mean synaptic gating variable (Nicola and Campbell, 2013). 
Mean-field systems enable one to conduct network level bifurca- 
tion analysis and hence to test different hypotheses about large 
network behavior. For example, Hemond et al. (2008) found that 
when uncoupled, the majority of hippocampal pyramidal neu- 
rons in region CA3 do not display bursting. This is contradictory 
to the observation that burst firing is ubiquitous in this region 
(Andersen et al., 2006, section 5.3.5). A possible explanation for 
these contradictory observations is that bursting is a network 
level phenomenon. This hypothesis has been tested using bifur- 
cation analysis of the mean-field system derived from a network 
of Izhikevich neurons (Dur-e-Ahmad et al., 2012). In particular, 
it was shown that for a network of identical all-to-all coupled 
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neurons fit to experimental data from Hemond et al. (2008) 
for CA3 pyramidal cells, bursting occurs for a large range of 
synaptic conductances and applied currents if the spike frequency 
adaptation in the neurons is sufficiently strong. 

Hemond et al. (2008) also observed that the pyramidal neu- 
rons in their study were heterogeneous, in particular, the neurons 
had different degrees of spike frequency adaptation. When the 
study of Dur-e-Ahmad et al. (2012) was extended to a network of 
two homogeneous subpopulations of hippocampal neurons with 
different degrees of spike frequency adaptation, the mean-field 
equations predict, and numerical simulations confirm, that the 
region in the parameter space where bursting occurs decreases in 
size (Nicola and Campbell, 2013). This would seem to indicate 
that adaptation induced bursting may not be robust to hetero- 
geneity, however, it is unknown how robust this network level 
bursting is to heterogeneity in different parameters. An extension 
of the mean-field system in Nicola and Campbell (2013) to large 
networks of heterogeneous neurons, is needed to fully analyze the 
robustness of adaptation induced bursting. 

The application of mean-field analysis to large networks of het- 
erogeneous neurons has been primarily limited to networks of 
one dimensional integrate and fire neurons with heterogeneity in 
the applied current (Hansel and Mato, 2001, 2003; Vladimirski 
et al, 2008). Hansel and Mato (2001, 2003) analyze a network 
of all-to-all coupled quadratic integrate and fire neurons con- 
sisting of two subpopulations: one excitatory and one inhibitory. 
They showed analytically that the tonic firing asynchronous state 
can lose stability either to a synchronous tonic firing state or a 
bursting state. Vladimirski et al. (2008) analyze a network of all- 
to-all coupled linear integrate and fire neurons subject to synaptic 
depression. This model was used to study network induced burst- 
ing in the developing chick spinal cord. Their derivation of the 
mean-field model is based on temporal averaging of the fast 
synaptic gating variable, in addition to the usual network aver- 
aging. This results in a model that only involves the distribution 
of slow synaptic depression variable and the distribution of fir- 
ing rates. Vladimirski et al. note that, for their model, increased 
heterogeneity tends to make population induced bursting more 
robust. They also note that one cannot understand the behavior 
of their network with a single slow, network averaged synaptic 
depression variable. 

More recently, Hermann and Touboul (2012) considered net- 
works with heterogeneity in the synaptic conductances. They 
compare heterogeneity induced by a distribution of parameters 
and heterogeneity induced by noise. For a firing rate (Hopfield- 
type) model, they derive a mean field representation of the 
situation with noise, which is a system of two ODEs for the 
mean and variance of the network average voltage. They show 
that increasing the strength of the noise, (which corresponds to 
the variance of the heterogeneity) causes a transition from qui- 
escence to oscillations in both the mean-field model and the full 
network simulations both with noise and a distribution of param- 
eters. Based on these results, they suggest a mean field model for 
a network of excitable Fitzhugh Nagumo neurons, which con- 
sists of coupled ODEs for the network mean voltage and network 
mean recovery variable. In both the mean-field model and full 
network simulations they observe a transition from quiescence to 



periodicity and then to chaos as the variance of the heterogeneity 
is increased. 

The motivation for the present paper is to explore the effect of 
heterogeneity in parameters on network induced bursting when 
adaptation is the primary negative feedback acting on the individ- 
ual firing rates. To this end, we introduce a set of mean-field equa- 
tions for networks of heterogeneous two-dimensional integrate 
and fire neurons. While the specific neural model we consider is 
for neurons with adaptation, our derivation is quite general and 
could be applied to other integrate and fire models. In contrast 
with (Vladimirski et al., 2008) our derivation of the mean-field 
model does not use temporal averaging thus we end up with a 
mean-field system which involves the network averaged synaptic 
gating variables as well as the distribution of adaptation variables. 
We also allow for heterogeneity in more than one parameter. Our 
approach is a generalization of that used for homogeneous net- 
works of two dimensional integrate and fire neurons (Nicola and 
Campbell, 2013), however, in the heterogeneous case it turns out 
that there are actually multiple mean-field models, as different 
assumptions can be made during the derivation. This leads us 
to three distinct mean-field systems, each derived under differ- 
ent assumptions, and used for different purposes. Together, these 
sets of equations allow us to do bifurcation analysis on large net- 
works, as in the homogeneous case. We show that the bifurcation 
structure of the heterogeneous network differs both qualitatively 
and quantitatively from the homogeneous network case. We dis- 
cuss the implications of this for network induced bursting in the 
hippocampus. 

When considering a homogeneous network, the mean-field 
variables are a good approximation for the variables of every 
neuron. However, this is not the case for a heterogeneous net- 
work. If the heterogeneity is large, then the neuron variables may 
be also widely distributed, rendering information about the first 
moments less useful. This also implies that the behavior of any 
individual neuron is less predictable with a mean-field system 
than it was in the homogeneous case. One of our mean-field 
systems addresses this problem, giving information about the 
distributions of the variables instead of just the mean. 

When considering a model for a specific heterogeneous net- 
work of neurons, one stumbling block is determining the dis- 
tribution of parameters. Estimates of the distribution can be 
made through direct intracellular recording of a sufficient num- 
ber of neurons and conventional measurements of the biophys- 
ical properties (membrane capacitance, voltage threshold, etc.). 
Unfortunately, this is a very time consuming and intensive pro- 
cess. What is needed is a way of measuring the biophysical param- 
eters of multiple neurons simultaneously. There are a few ways 
to sample multiple neurons such as multi-unit recordings using 
tetrode arrays, or two-photon microscopy techniques (Buzsaki, 
2004; Crewe et al, 2010). However, these techniques typically 
only tell us about spike times of large (dozens to hundreds of 
neurons) networks (Buzsaki, 2004). While an impressive accom- 
plishment, this still does not tell us anything directly about the 
biophysical properties of the neurons that caused those spikes. In 
this paper we use a mean-field system to determine an approx- 
imate distribution of firing rates for a network given a known 
distribution of parameters. This is an unusual state of affairs for a 
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mean-field system, as these kinds of systems seldom give infor- 
mation about entire distributions. More importantly, however, 
using a mean-field system, we can invert a distribution of steady 
state firing rates (which can be obtained from multi-unit record- 
ings) to obtain a distribution of parameters. In fact, this can be 
done at the individual neuron level, to determine the parameter 
value for any particular neuron. This allows one to estimate dif- 
ferent biophysical parameters, which are difficult to measure at 
the network level, using easy to measure firing rate distributions. 
However, the assumptions required for the numerical accuracy of 
the estimation are fairly strong. 

The plan for our article is as follows. Section 1.1 introduces 
the general class of adapting two-dimensional integrate and fire 
neurons used in our network models. This class was introduced 
by Touboul (2008), who also completed the bifurcation analysis 
of the single, uncoupled neuron. Population density methods are 
briefly introduced in section 1.2, as a population density equation 
serves as an intermediate step to obtain our mean-field models. 
Section 2 begins with a review of mean-field theory and the equa- 
tions for the homogeneous network. This is followed, in sections 
2.1-2.3, by the derivation of the three mean-field systems for the 
heterogeneous network. A comparison of numerical simulations 
of these mean-field systems and the full network is the subject 
of section 2.4. Applications of mean-field theory to networks 
with a single heterogeneous parameter can be found in section 
3 including bifurcation analysis (section 3.1), distributions of 
parameters and firing rates (section 3.2) and using mean-field 
theory for parameter estimation from firing rate data (section 
3.3). Applications of mean-field theory to networks with multiple 
sources of heterogeneity are included in section 4. A discussion of 
our work and its implications can be found in section 5. 

2. MATERIALS AND METHODS 

2.1. NON-LINEAR INTEGRATE AND FIRE NEURONS 

We consider a network of two-dimensional integrate and fire 
models of the form 



v = F{v) — w + I 
w = a(bv — w), 



(1) 
(2) 



where v represents the non-dimensionalized membrane poten- 
tial and w serves as an adaptation variable. Time has also been 
non-dimensionalized. The dynamical equations (1, 2) are supple- 
mented by the following discontinuities 



(^spike) 



V l f spike) 
W (/spike ) 



^reset > 



(fspike ^) 



+ W\ 



(3) 



jump ■ 



This particular notation was formally introduced by Touboul 
(2008), along with a full bifurcation analysis of this general family 
of adapting integrate and fire neurons. Members of this fam- 
ily include the Izhikevich model (Izhikevich, 2003), the adaptive 
exponential (AdEx) model (Brette and Gerstner, 2005; Naud et 
al., 2008) and Touboul's own quartic model (Touboul, 2008). 

The methods of this paper can be applied to a network of any 
particular neuron belonging to this general family, and thus all 



derivations are done for this model. For the numerical examples, 
however, we only consider the Izhikevich neuron. In dimensional 
form this model is: 



CVi = k(V, - V T ){Vi - V R ) - Wj + 7 apP , 
TKVj - Vr) - Wj 
tw 



W; 



spike J 



W ! feke) = W .( f spike) 



+ Wj 



(4) 
(5) 

(6) 



jump, z> 



In dimensionless form, this model is given by Equations (1-3) 
with F(v) = v(v — a) in addition to dimensionless versions of 
the discontinuities (Equations 6). We will use uppercase letters 
for dimensional variables and lower case for their dimensionless 
counterparts. The application to other neural models is straight 
forward, see Nicola and Campbell (2013) where the homoge- 
neous mean-field theory has been derived and tested for both the 
AdEx and the Izhikevich models. 

Networks of these neurons can be coupled together through 
changes in the synaptic conductance. The synaptic conduc- 
tance of post-synaptic neuron i due to presynaptic neurons j = 
1, 2 -NTs given by 



?i(t) = giSi(t) 



N 



(7) 



where g, denotes the maximal synaptic conductance of neuron i 
and s,(f) denotes the total proportion of postsynaptic ion chan- 
nels open in the membrane of neuron i. The time dependent 
variable s«(f) represents the proportion of postsynaptic ion chan- 
nels open in the membrane of neuron i as a result of the firing in 
neuron j. 

The changes in Sy(t) that occur after a spike are often mod- 
eled as transient pulses. For example, if neuron j fires its fcth 
action potential at time t = t; k, then the variable sy(t) at time 
t is given by 

tj,k<t 

There are different functions proposed for E(t) in the literature 
including the simple exponential, the alpha synapse and the dou- 
ble exponential. We primarily consider the simple exponential 
synapse 



E(t) 



J jump 



exp 



which is governed by the ordinary differential equation 
dsdt) 



dt 



+ Sjump 2l 8 (f - t h k ) 



(9) 



(10) 



In the rest of the paper, we assume all-to-all connectivity and 
that the synaptic parameters Sj um p and x s are the same for every 
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synapse. In this case we may set s;(f) = s(f) for all i, as each post- 
synaptic neuron receives the same summed input from all the 
presynaptic neurons. Then, using Equations (7) and (10), the net- 
work of all-to-all coupled neurons that we consider is given by the 
following system of discontinuous ODE's: 

n = F(v ; ) - w, + I, + g,s(t) (E r - v.) , (11) 
w; = a\{bv{ — wi), (12) 

* = -f + s i? EE *<*-**>. d3) 



j=\tj_ k <t 



v <( f >, It) = Vak 



Vi(t^jfe) = Preset, 



(14) 



jump • 



forj = 1,2, ...N. 

In the examples, we consider one or more parameters as the 
sources of heterogeneity. However, to simplify the notation in the 
derivations, we use the vector P to represent all the heterogeneous 
parameters. Then, denoting the state variables v and w as the vec- 
tor x, we can write the equations for the individual oscillator as 



x = G(x, P, s) 



ch(x,M 

G 2 (x, P) 



(15) 



Given a specific heterogeneous parameter, Gi and G 2 may not 
depend on P, or all of the components of p. However, for the sake 
of simplicity, we include the dependence in both equations. 

Our numerical examples are restricted to the Izhikevich neu- 
ral model, and we primarily consider the driving current J, of 
each neuron, the synaptic conductance g; and the adaptation 
jump size, Wj ump ,,- as the source of heterogeneity. However, the 
mean-field equations we derive can be applied to any of the 
two-dimensional adapting integrate and fire models, with any 
heterogeneous parameter or set of parameters. 

Finally, we note that in many applications b is a small param- 
eter, and thus the bv term can be dropped in G 2 - We do this in 
all our numerical studies. However, one can still derive appropri- 
ate mean-field equations if this term is present (see discussion in 
Nicola and Campbell, 2013), and thus we have left the term in the 
derivations. 

2.2. THE POPULATION DENSITY EQUATION 

The population density function, p (x, t) determines the density 
of neurons at a point in phase space, x, at time t . Consider first the 
case of a homogeneous network, i.e., all the oscillators have the 
same parameter values, denoted by p. In the limit as N —> 00, one 
can derive the following evolution equation for the population 
density function: 



3/0 (x, f) 
dt 



-V • J(x, P, 5, f) 



where J is given by 

Kx, P, s, t) = G(x, p, s)p{x, t) = (f, J w ) . 



(16) 



(17) 



and must satisfy the boundary condition 

7 V (Vpeak, W, P, S, f) = J V (v [eset ,W+W) ump , P, S, t). (18) 

In the same limit, the differential equation for s converges to 



+ / />peak,W,S,P, t)dw (19) 

Is JW 

where the integral term is actually the network averaged firing 
rate, which we denote as (R(t)}. Derivations of Equation (16) 
can be found in various sources (Nykamp and Tranchina, 2000; 
Omurtag et al., 2000). 

Equation (16) is frequently referred to as the continuity equa- 
tion and it has various applications besides its use as an inter- 
mediate step in mean-field reductions. For example, the equation 
has been used to determine the stability of the asynchronous 
tonic firing state by various authors (Strogatz and Mirollo, 1991; 
Abbott and van Vreeswijk, 1993; van Vreeswijk et al., 1994; 
van Vreeswijk, 1996; Hansel and Mato, 2003; Sirovich et al., 
2006). These papers predominantly consider homogeneous net- 
works of linear integrate and fire neurons. The exception is the 
work of Hansel and Mato (2003) which considers heterogeneity 
in the applied current. One can study stability of various fir- 
ing states using spectral analysis or other analytical treatments 
of this equation (Strogatz and Mirollo, 1991; Abbott and van 
Vreeswijk, 1993; van Vreeswijk, 1996; Knight, 2000; Sirovich 
et al, 2000, 2006; Hansel and Mato, 2003). However, these 
approaches are too complicated for the models we consider in 
this paper. 

Now consider a heterogeneous network where the param- 
eters vary from oscillator to oscillator, but are static in 
time. Then one can rewrite the equations for the individual 
oscillator as 



v, = Gi(x„ P„ s), 
Wj = G 2 (Xi, P ; ), 

Pi = 0. 



(20) 
(21) 
(22) 



In this case the flux contribution due to P is 0 and the evolution 
equation for the network is given by 



dp(x, p, t) 
dt 



-V • J(x, p, s, f) 



(23) 



The density now has the vector of parameters, P, as an indepen- 
dent variable. The flux consists of the vector (J v , J w , 0), with P 
as an independent variable, as opposed to a fixed constant. If the 
parameters are time varying however, the final component of the 
flux will be non-zero. The equation for s is also different in the 
heterogeneous case: 



+ Sjump [ [ /Vpeak, W, 5, P', t) dw d$ . (24) 

Jw Jr 



While the evolution equation (23) is an exact representation 
for the network in the large network limit, it is difficult to 
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work with analytically. Additionally, as the dimensions of the 
PDE become large, it becomes difficult to find numerical solu- 
tions efficiently (Ly and Tranchina, 2007). However, mean-field 
reductions of the network can be used to reduce the popula- 
tion density PDE to a system of non-linear switching ODEs that 
governs the moments of the distribution. Unlike the PDE, the 
system of ODEs is tractable using bifurcation theory, at least 
numerically. Furthermore, we show that in the heterogeneous 
case, the resulting mean-field systems can yield more informa- 
tion than just the type of bifurcation that the network can 
undergo. 

2.3. MEAN-FIELD THEORY 

In the homogeneous case, the mean-field system of equations for 
an all-to-all coupled Izhikevich network was derived in Nicola 
and Campbell (2013). We present here a quick summary of 
this derivation. In order to derive a mean-field system of equa- 
tions, one first needs to reduce the PDE for p (x, t) = p (v, w, t) 
by a dimension. This is done by first writing the density in its 
conditional form: 



p{v, w, t) = pv(v, t)pw(w\v, t) 



(25) 



and then integrating the continuity equation with respect to w. 
This yields the one dimensional PDE 



dp v (v, t) 



9Gi(v, (w\v), s)pv(v, t) 



dj(v, (w\v),s, t) 



dt dv 
where the flux has been redefined to 



dv 



J(v, (w\v),s, t) 



Jw 



J (v, w, s, t) dw. 



(26) 



(27) 



One can now make a first order moment closure assumption, 
(w\v) = (w), and derive an approximate ODE for (w), which 
yields the system 



9 3 
p(v, t) = - — ((F(v) - (w) +I + g(E r - v)s) (p(v, f))) 



dt 



dv 

. b(v) - (w) 

( W ) = 1- Wj ump /(v peak , (w),s, t) 



S = 1" 5j um p7(v peak , (w},S, f) 



where the subscript on the density function has been dropped for 
convenience. The details and validity of the first order moment 
closure assumption that is used can be found in Ly and Tranchina 
(2007). We note, however, that the work in Ly and Tranchina 
(2007) was primarily with leaky integrate and fire neurons, as 
opposed to the two-dimensional adapting class we consider here. 
However, it is a necessary assumption to proceed analytically. If 
we assume that the adaptation time constant, x w = - is large, one 
can apply a quasi-steady state approximation to derive a system of 



switching ODE's for (w) and s: 
b(v) - (w) 



+ W\ 



lump 



(28) 



(R t (t)) 



s 

[/- 



'jump 



(RM) 



(29) 



dv 



V F(v)-{w) + I + g(e r - 



v)s\ 



■■ H((w),s) > 0 
: H((w),s) < 0 



(30) 



The switching manifold for the system, H((w), s) is given by: 
H«w), s) = 1- (w) + min(F(v) + g(e r - v)s). 



(31) 



Note that H((w), s) depends on the parameter(s) of the model, 
and thus for the heterogeneous case, we make this dependence 
explicit by writing H((w), s, |3). As the computation for (v) is 
somewhat lengthy and is only outlined in the discussion of Nicola 
and Campbell (2013), we have placed it in Appendix A. Note 
that this approach is similar to temporal averaging of a fast volt- 
age equation assuming slow synaptic and adaptation currents, 
as outlined in Ermentrout (1998) and Ermentrout and Terman 
(2010). 

For the Izhikevich neuron, Equations (30, 31) become 



(R,(t)) 



. It* 



dv 



V v(v—a)—(w)+I + g(e r — v 



J" 1;H(< 



w), s) > 0 



H({w),s) < 0 



(32) 



H({w),s) = I-(w) 



(a + gs) 2 



+ 1 



(33) 



Note that in this case, we can evaluate (Rj(t)} explicitly: 

0 : H({w), 5) < 0 



v peak 2 
V'- 1 *((«■>, s) 



in addition to an approximation to ■ 



<v) = 



WO) 



log 



(v p eA-'^ £ ) 2 +H((M>),5) 

(Vreset-^) 2 +H((w),s) 

V-H((w>,s) 



H((w), s) > 0 



: H((w), s) < 0 
(34) 

A comparison of solutions of these equations and the full network 
are shown in Figure 1 for both the tonic firing and the burst- 
ing case. Note that during tonic firing, the steady state firing of 
the individual neurons appears to be asynchronous, while during 
bursts the neurons separate into two synchronous subpopula- 
tions that fire out of phase with one another (see Figures 1E,F). 
While the mean-field system is more accurate when the neurons 
all fire asynchronously, the synchronization in the bursting state 
does not appear to be a substantial source of error in determin- 
ing the mean-values and the resulting qualitative behaviors. The 
asynchronous firing in the non-bursting region is consistent with 
previous work on the stability of asynchronous states with exci- 
tatory coupled class one oscillators (Abbott and van Vreeswijk, 
1993). 
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FIGURE 1 | Numerical simulation of a homogeneous network of 1000 
Izhikevich neurons, with parameters as follows. (A,C,E) l app = 4500 
pA, g syn = 200 nS (B,D,F) / app = 3500 pA, g syn = 200 nS. The rest of the 
parameters can be found in Table 1. Simulations of the mean-field 
equations (in red) and the mean values of the corresponding full network 
simulations (in blue) showing (A) tonic firing and (B) bursting. In this 
and the following figures (W(t)) is the network mean adaptation variable 
in dimensional form and (g(t)) = g S ynS(f) is the network mean synaptic 



conductance. (B,D) Raster plots for 40 randomly selected neurons from 
the network simulation in (A,C). (E,F) are the last 100 ms of the raster 
plots in (CD), respectively. The mean-field equations are fairly accurate 
both when the network is tonically firing and when it is bursting. For 
the tonic firing case the neurons fire asynchronously at steady state, 
while in the bursting case, the neurons seem to align into two 
synchronously firing subpopulations firing out of phase with one another 
during the burst. 



This system of equations is valid when x w 3> O(l), however, 
the magnitude of x 5 is also significant. While in the original 
derivation of Nicola and Campbell (2013), x s = 0(x w ) was sug- 
gested as a criterion for validity, this is not actually necessary. One 
merely requires that x 5 not be significantly smaller than O(l), the 
time scale of the PDE. The reason for this is that if the time scale 
of the ODE for 5 is smaller than that of the PDE then the quasi- 
steady state approximation must be applied to the ODE for s as 
well. The requirements on the time constants are carried forward 
in the heterogeneous case. 



In our models, the timescale of the ODE for s is typically 
between that of the PDE and that of the ODE for w, thus we 
have not applied the quasi-steady approximation to s. Applying a 
quasi-steady state approximation to both s and the reduced PDE 
yields a more compact system, which is just an ODE for ( w) , how- 
ever, the analysis does not get any simpler. The reason for this is 
two-fold: the ODE for (w) remains non-smooth and the firing 
rate now has to be implicitly solved at each time step. Thus, it is 
more convenient to apply the quasi-steady state approximation 
only to the partial differential equation. 
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When parameter heterogeneity is added into the mix, it turns 
out that there are multiple "mean-field" systems of equations that 
can be derived, by applying different assumptions on the condi- 
tional moments. We outline three different assumptions that can 
be made and derive the resulting system of mean-field equations 
in each case. 

2.3. 1. Mean-field system I 

We begin by writing out the density function in the conditional 
form 

p(x, P, t) = p x (x, 0ppO[*. 0 OS) 
The continuity equation is then given by 
3(px(x, f)/°p(PI*, 0) 



dt 



v.j(*, s,M). 



(36) 



Simple integration with respect to ft yields the reduced continuity 
equation 

dp x (x, t) 



dt 



-V-J(x, s, (p|x),t). 



(37) 



This step is valid for all the non-dimensionalized models we 
consider as they are all linear in their dimensionless parame- 
ters (see Touboul, 2008). The flux has also been redefined upon 
integration to 

(J V J W ) = p s {x, t) (Gi(*, (0|x), s), G 2 (x, (P|x))) . 

We now apply the moment closure assumption (P|x) = (P) to 
yield the following PDE: 



dp x (x, t) 
dt 



-V-/(x,s, <P),t). 



(38) 



It should be clear that this is equivalent to the continuity equation 
for a homogeneous network with parameter values fixed at (P). 
Thus, the associated mean-field system is identical to the homo- 
geneous case, only with the parameters fixed at (P). This is the 
simplest assumption one can make in the heterogeneous case. For 
example, if we treat J as the source of heterogeneity for a network 
of Izhikevich neurons, with distribution pi(I), then the resulting 
mean-field system is 



(w)' = 



b(V) - (W) 



ARM) 



+ Sjump (#;({)) 



<«;«> = 



(A 



dv 



V v(v-u)-(w) + (I)+g(e, 
0 



H((w),s, (/» = (I)~(w)- 



(a + gs) 1 



(39) 
(40) 

H((w),s, (J))>0 
(41) 

H((w),s, (/» <0 
(42) 



(v) 



(M>) 

2 



Y+H((W),S, (/)) 



\(Vw«-^) z +H((w),s, (/)) 

:H((w),s, (7))>0 

^-V-H((w>, 5, (/» 
:H((w>, s, (J)) < 0 



(43) 



Note that 7 in Equations (32, 33) has been replaced by (I) in 
Equations (41-43). We treat this system as the baseline mean-field 
model for comparison purposes, in addition to direct numerical 
simulations of the network. We denote this system of equations as 
mean-field one (MFI). We should expect this system to be an ade- 
quate approximation to the actual network for narrowly centered 
distributions of the parameter heterogeneity (small values of the 
variance, op). 

This set of differential equations is representative of a common 
approach taken when fitting actual neurons. In this approach, 
multiple estimates of parameters or measurements taken from 
multiple neurons are averaged to yield a single parameter value, 
which is really the mean parameter value, (P). Simulations of 
homogeneous, large networks are then run with the parameters 
fixed at their mean values. As we shall see in subsequent sections, 
the behavior of a simulated heterogeneous network can differ 
substantially from that of MFI. 

2.3.2. Mean-field system II 

To derived our second mean-field system, we begin by writing the 
density function in the alternative conditional form 

p(v, w, P, f) = Pw (w, f|P, v)p v (v, t|P)pp(P). (44) 

Next we integrate the continuity equation with respect to w. This 
yields the following system 



— Pfs(P) = - / 

dt Jw \ 



s, P, t) 9/ w (v, w, s, p, t) 
dw 



dw 



3 

dv 
3 

3v 



7(v, (w|v, p>,s,p,f)-/ w (v, w, s,P,f)|aw 
J(v, (w\v, P>, s, P, f), (45) 



where the last term vanishes as J w is assumed to be vanishing on 
the boundary, and 



/O, (w|v,P),s, P, t) = / J v (v,w,s, P, t)dw. (46) 
Jw 

We now make the first order moment closure assumption 
(w\v, P) = (w). Then to complete the system, we must derive a 
differential equation for (w) : 



! dwdv 



dp(v, w, P, f) 

> 

dt 

ii' I 1 ) d^ dwdv 

dw dv 



Jv Jw J$ 

--II L- 

Jv Jw Jp 

= / / G 2 (v,w,$)p(v,w,$,t)d$dwdv 
Jv Jw Jp 

- / / w(/ v (v peak , w, s, p, f) - / v (v reset , w, s, p, t)) dp dw 
Jw Jp 

= (G 2 (v, w, p)> - [ [ w(/ v (v peak , w, 5, p, f) 
Jw Jr 
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-/^(Vpeak, W - Wj ump , 5, P, f)) 



(47) each neuron has a total input current given by J; + 7 syn , where 
J syn is the synaptic current given by the network coupling. It 



(G 2 (v, w, P)) + / / Wjump/Vpeak, W ' S ' P> f ) dw d P foUows that the network averaged firing rate is approximately a 
JRJW function of / syn 



+ 0(<mp) 



G 2 ((v>, (w), (p» + [ w jum p/(v peak , (w), s, p, f) dp. (48) 
Jr 



Note that we have made the approximation (G2(v, w, P)) = 
G2((v), (w), (P)) in addition to dropping the 0(w? ump ) terms. 
Additionally, the substitution in Equation (47) comes from the 
boundary condition (Equation 18). 

Applying a quasi-steady state approximation to the PDE 
(Equation 45) yields the following equation for the steady state 
voltage independent flux, J({w), s, P): 



7«w),s, 



}o if H((w),s, P) < 0 

(49) 

We interpret the ratio J((w), s, P)/pp(P) as the parameter depen- 
dent (or conditional) network averaged firing rate, (_R,(£)|P), 
based on the fact that 



J((w),s, P) dp« (R,(t)). 



In other words, the distribution of parameters induces a distribu- 
tion of firing rates across the network, and the network averaged 
firing rate is the mean of the distribution. 

In summary, the resulting mean-field equations are given by: 



(w)' 



b(v) - (w) 



+ f w jump (^(0IP>Pfi(P)dp(50) 
Jr 



'jump 



/W)IP>Pf5(P)dp (51) 

Jr 



(Ri(f)IP) 



[fvW^y 1 :H((w>,5,P)>0 (52) 
0 :H((w),s, P) < 0 

H((w) , s, P) = I — (w) + min(F(v) + g(e r - v)s) (53) 



(v) 



/ (v|P)pp(P) 

Jr 



rfp 



(54) 



where the forms of Gi(v, (w), s, P) and H((w), s, P) depend on 
which specific neural model is used, and the equation for (v|P) 
can be found in Appendix A. Note that the distribution of fir- 
ing rates is not computed explicitly in these equations, only the 
conditional firing rates, (R;(f)|P), are computed. However, we 
show in section 3.2.2 that a distribution for the steady state firing 
rates of the network can be computed using (R;(f)|P). We refer to 
Equations (50-54) as mean-field two (MFII). 

It appears that MFII adds some smoothness to the non-smooth 
MFI equations. This can be easily seen by taking a heteroge- 
neous background current to each neuron, In this situation, 



If 7 syn is treated as parameter, this equation can be evaluated with 
differing standard deviations for the normally distributed input 
current. When this is done we see that the F(I) curve becomes 
smoothed out as the standard deviation increases. This is shown 
in Figure 2. 

Additionally, MFI and MFII also differ in the order in which 
the integrations are carried out. In MFI, we integrate with 
respect to P first, and then apply the first order moment clo- 
sure assumptions (P|x) = (P) and (w\v) = (w). In MFII, we 
integrate with respect to w first, and then apply the moment clo- 
sure assumption (w\v, P) = (w). Furthermore, if (_R,(f)IP) does 
not actually depend on the heterogeneous parameter P, such 
as when the heterogeneity is in Wj U mp) then MFI and MFII are 
identical. 

The first order moment closure assumption used here can 
be weakened. This leads to the "mean-field" system in the next 
subsection, which is a different kind of system than MFI and 
MFII. 

2.3.3. Mean-field system III 

Suppose that instead of assuming that (w\v, P) = (w), we make 
the weaker assumption that (w\v, P) = (w|P). It turns out that 
this assumption yields a PDE, even when one makes the quasi- 
steady state approximation, as we now show. Applying this weaker 





= 0 pA 




= 50 pA 




= 100 pA 


— 0/ 


= 250 pA 




0 500 1000 1500 2000 

Current (pA) 

FIGURE 2 | The firing rate curve HI) for a Gaussian distributed 
background current in the (R,it)) response curve for the homogeneous 
network. The F(l) curves are plotted for increasing values of o/ which 
smooths out the square root type non-smoothness at the onset of firing. 
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moment closure assumption to Equations (45) yields the follow- 
ing simplification of the continuity equation: 



3pv(v, f|p) 



Pp(P) 



-7(v, <w]p),s,p,f). 



dt 9v 
Application of the quasi-steady state approximation now yields 



(55) 



/(V,MP),5, p) 



0 :H((w|P>,s, P) < o' 



H((HP>, 5, P) = 7 - (HP) + min(F(v) + g(e r - v)s). 

An equation for the time variation of (HP) can be derived in a 
similar manner to the last section, yielding the following mean- 
field system: 



9(HP> fo(v|P)-(w|P) 



dt 



ds 
dt 

(WW) 



+ S: 



'I Limp 



+ H>™p(- R >«IP> 

f (K,(f)IP>Pfs(P)dp 

JR 



:H«w|P),5, P) <0 



(56) 
(57) 

(58) 



Note that (w) can be computed via: 



(w) 



[ (HP)/Op(P) dP- 
Je, 



(59) 



The equation for (v|P) can be found in Appendix A. We denote 
this system as mean-field three (MFIII). Note that the equation 
for (HP) is actually a PDE. This is due to the fact that the con- 
ditional moments, (HP), (-RIP) an d MP) are functions of both 
time and the "spatial" variable p. This partial differential equation 
is easier to deal with than most PDEs as it has no spatial deriva- 
tives, however, the right hand side of the differential equation is 
non-smooth. 

While this system should be more accurate than mean-field II, 
it has the drawback of being more difficult to analyze. The depen- 
dence on P forces one to discretize over a mesh in P in order 
to work numerically with this system. This approach, typically 
referred to as the method of lines in the literature, is often used 
to solve PDE's with no spatial derivatives. We use this approach to 
numerically simulate Equation (56). In order to compute the inte- 
grals in Equations (57-59) using the method of lines, we choose 
a grid that is non-uniform and generated with the density func- 
tion pp (P) . The integrals are subsequently replaced with averaging 
over the entire grid, which is precisely a Monte-Carlo method for 
estimating the integrals. 

Numerical bifurcation analysis of MFIII is more difficult as 
it is a PDE. In principle it is possible to do numerical bifurca- 
tion analysis on a PDE by discretizing and analyzing the result- 
ing large system of coupled ODE's (Ko and Ermentrout, 2007). 
However, when we tried this approach on MFIII it proved to 
be too numerically intensive, and required a lengthy period of 



time for convergence of the numerical methods used for con- 
tinuation. Additionally, the system of ODE's is still non-smooth, 
which causes problems with most numerical continuation soft- 
ware. However, as we shall show later, an approach that yields 
similar information to direct bifurcation analysis can be used with 
MFIII. 

3. RESULTS 

3.1. NUMERICAL SIMULATIONS 

We carried out numerical simulations of the full network model 
with a large number of neurons and the corresponding mean 
field systems. The large network simulations are carried out on 
the dimensional version of the equations. The results are pre- 
sented in terms of the network mean adaptation, (W(t)}, which 
is the dimensional version of (W), the network mean synaptic 
conductance, (g{t)) = gsynSit), and the dimensional parameters 
described in Table 1. Since the mean field systems are given in 
terms of the dimensionless variables and parameters, results from 
mean field simulations are converted to dimensional form for 
comparison with the full network simulations. 



Table 1 | The values of the model parameters and variances of the 
distributions used in this paper. 



Dimensional parameters 

C 250pF 
k 2.5nS/mV 
V R -65 mV 

V T Vr + 40 - r 



= 41.7mV 
30 mV 
-55 mV 
200 pA 

200 ms 
-1 nS 

1000-5000 pA 
0-600 nS 



Dimensionless parameters 



''peak 

Weset 

^ump 

XW 

'1 



a = 1 + 
a = 1 + 

^peak — 



\Vr\ 
\V R \ 



1 + 



^peak 
Weset 



app 



£?syn 



N 
oi 
"g 

°d 

m (mixing 
parameter) 



^ w k\V R \ 



k\V R \ 



4 ms 

0.8 

1000 

0-1000 pA 
50 nS 
50 pA 
0-1 



k\V H \ 2 

0syn 

' Wr\ 

_ tsyn^l^fl 

c 



0.6215 

0.6215 

1.461 

0.1538 

0.0189 

0.0077 

-0.0062 

0.0776-0.3333 

0-3.6923 

2.6 



These values apply unless otherwise indicated. Rheobase for the dimensional 
parameter values is l rh = WOOpA. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



December 2013 | Volume 7 | Article 184 | 9 



Nicola and Campbell 



Mean-field models for heterogeneous networks 



Recall that simulations of a homogeneous network and the 
corresponding mean-field system are shown in Figure 1. Note 
that the network undergoes a bifurcation from tonically firing to 
bursting as the amount of applied current / app is decreased, with 
all other parameter values held fixed. Further simulations show 
that if 7 app is decreased below 7 r h then all neurons in the network 
are quiescent (non-firing). 

To determine and compare the validity of the three mean- 
field systems we derived for the heterogeneous networks, we 
have run a series of numerical simulations of these systems and 
of an actual network containing 1000 neurons. The parame- 
ter values for the individual neurons can be found in Table 1. 
They are based on those given in Dur-e-Ahmad et al. (2012) 
which were fit to data for hippocampal CA3 pyramidal neu- 
rons from Hemond et al. (2008). These are the parame- 
ter values we use for the rest of this paper, unless otherwise 
indicated. 

As a starting point, we consider heterogeneity only in the 
applied current. The distributions are assumed to be normal 
with mean (7) and standard deviation 0j. We varied the values 
of the mean and standard deviation and found that the accuracy 
of the mean-field approximations depends on where the mean is 
relative to the different bifurcation regions and on the size of the 
standard deviation. 



A 6000 




0 200 400 600 800 1000 

Time (ms) 



FIGURE 3 | Numerical simulations of a network of 1000 Izhikevich 
neurons with parameters as in Table 1, except g syn = 200 and the 
applied current which is normally distributed with mean and variance as 
follows. (A) (/ app ) = 5000 pA, a, = 2000 pA. (B) </ app ) = 3000 pA, 
a I = 200 pA. (C) (/ app ) = 3000 pA, ct, = 500 pA. (D) (/ app ) = 3000 pA, 
o/ = 2000 pA. Blue is the network average of a given variable, red is MFI, 



As for the homogeneous network, the heterogeneous net- 
work undergoes a bifurcation from tonic firing to bursting as the 
amount of current applied to the individual neurons is decreased, 
with all other parameters held fixed. This can be seen in Figure 3 
where the bifurcation with decreasing (7 a pp) is shown. As (7 a pp) is 
decreased below 7 r h, there is a bifurcation to quiescence. We will 
not discuss this latter bifurcation in detail, as we are primarily 
interested in analyzing the transition from tonic firing to bursting. 

Note that the bifurcations described above only occur in the 
mean sense. Since the current values are normally distributed, 
there is non-zero probability that some neurons receive large 
enough or small enough current to be in a state other than that 
corresponding to the value of (7 app ). For small enough standard 
deviations, very few neurons in an actual finite network are likely 
to have behavior different from the mean. However, for large stan- 
dard deviations, a sizable proportion may not follow the mean 
behavior. 

Given this knowledge of the different qualitative behaviors of 
the network, we can see how the mean-field systems compare. 
For tonic firing (Figure 3A), even when the standard deviation 
is large, the mean-field systems approximate the network means 
(g( f )) an d ( W(f)) very well. However, when the network is burst- 
ing, with (7 a p P ) > -T r h, we see a difference as to which mean-field 
system is superior. For small values of 07, we have numerically 
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green is MFII and black is MFIII. In this region, the mean-driving current is 
away from rheobase, (/ app ) 2> / rn - All three approximations are quantitatively 
and qualitatively similar for small to intermediate sized variances in the 
distribution of currents. For small variances, MFI is the most accurate and for 
larger variances, MFIII is the most accurate. For large variance, MFII 
bifurcates back to tonic firing earlier than MFI and MFIII, as seen in (D). 
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FIGURE 4 | Numerical simulations of a network of 1000 neurons with 
parameters as in Table 1, except g syn = 200 and the applied current 
which is normally distributed with mean and variance as follows. 

(A) (/ app ) = 1200 pA, cr, = 200pA. (B) (/ app ) = 1200 pA, o, = 500pA. (C) 
Capp) = 1200pA, cr, = 1000pA. (D) (4 PP ) = 1200 pA, cr, = 2000pA. Blue is 
the network average of a given variable, red is MFI, green is MFII, and 
black is MFIII. In these simulations, the mean-driving current is close to 



(and over) the rheobase. In all cases, MFI is the least accurate. This is 
because it depends only on (/ app ). When {/ app > = 0(/ r h), even for small 
variance, many of the neurons have / < l,h and may not spike at all. 
(A,B) For small values of 07, all three approximations are qualitatively and 
quantitatively accurate. (C,D) For larger variance, cr, = 0(/ r h), only MFIII is 
qualitatively and quantitatively accurate. In this case, MFII bifurcates early 
to tonic firing. 



found that mean-field I is superior to mean-field II and III, 
however all the systems are quantitatively and qualitatively accu- 
rate (see Figures 3B,C). However, for larger values of 07, the 
amplitude error of MFIII is the smallest, and MFII is the worst 
approximation as it bifurcates to tonic firing prematurely (see 
Figure 3D). 

When (-f a pp) is close to 7 r h, we see even stronger differences 
between the three mean-field systems. For small to intermedi- 
ate standard deviations, MFII and MFIII are clearly superior 
to MFI, having a smaller amplitude and frequency error (see 
Figures 4A,B). However, for larger values of cr; as shown in 
Figures 4C,D, only MFIII is a qualitatively and quantitatively 
accurate representation of the behavior of the network. The 
amplitude and frequency error of MFI are very large, and MFII 
again bifurcates prematurely to tonic firing. 

One should note that for (/ a pp) = 0(7 r h) and for large values 
of cr/, the network can undergo a period doubling bifurcation. 
This is shown in Figure 5. The large standard deviation in the 
current forces different neurons into different regimes, such as 
tonic firing, bursting, alternate burst firing and quiescence as 
seen in Figure 5A. During a burst the heterogeneity causes the 
neurons to fire asynchronously. Neurons with higher applied 
current fire followed by those with lower applied current. See 



Figure 5B. A small subpopulation of neurons with low applied 
current are alternate bursters (i.e., burst with twice the period of 
the rest of the bursting neurons). This appears as a period dou- 
bled limit cycle in the mean variables of the network, as seen 
in Figures 5C,D. Only MFIII is able to approximate the period- 
doubled limit cycle with any degree of accuracy, as shown in 
Figures 5C,D. Period doubling bifurcations are well known for 
their capability of inducing chaos. Given that MFIII accurately 
represents the period doubling bifurcation, it may be able to repli- 
cate any potential chaotic behavior. However, we leave further 
investigation of this interesting behavior for future work. 

To summarize, all the mean-field systems are valid for tonic fir- 
ing parameter regimes, and MFI is valid for all parameter regimes 
with small or, except for (/ a pp) = 0(J r h). Mean-Field II and III 
are valid for bursting with (-f app ) S> J r h> ana MFIII is the only 
valid approximation for {I app ) = 0(J r h). Thus, when (/ app ) is 
large we maybe able to use MFII to determine the type of bifurca- 
tion^) involved when a heterogeneous network transitions from 
tonic firing to bursting and the location in parameter space of 
the bifurcation curves. Note that when the mean network behav- 
ior undergoes a bifurcation from a tonic firing steady state to a 
bursting oscillation, this does not indicate that the entire net- 
work of neurons is bursting, or tonically firing. However, we will 
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FIGURE 5 | Period doubled limit cycle in the heterogeneous network and 

in MFIII. The network consists of 5000 neurons, with parameters as in 
Table 1, except g syn = 200 and the applied current which is normally 
distributed with mean (/ a pp) = 1100 pA and variance ct; = 2000 pA. (A) Raster 
plot of 100 randomly selected neurons of the network arranged in order of 
increasing current. Individual neuron behaviors include burst firing, alternate 



burst firing, tonic firing and quiescence. (B) Close-up of raster plot. The 
neurons appear to fire in "traveling waves" during a burst. (C) Comparison of 
the mean variables of the large network simulation and the simulations of the 
mean field systems. Only MFIII is able to reproduce the period doubling 
behavior. (D) Comparison of the "phase portrait" of period doubled limit cycle 
for MFIII and the mean variables of network. 



show how to use MFIII to determine what proportion of neurons 
display the different types of behavior, given a specific parameter 
regime and level of heterogeneity. 

In addition to simple heterogeneity using unimodal distribu- 
tions, one can also apply the same three mean-field equations to 
networks where multiple subpopulations exist. However, unlike 
previous attempts at modeling networks with multiple subpop- 
ulations, we do not generate discrete coupled subnetworks with 
different fixed values of the parameters in each subnetwork. 
Instead we use a smoother approach where the networks have 
distributions of parameters with multiples modes indicative of 
multiple subpopulations. This can be easily done through the 
processing of mixing unimodal distributions (see Appendix C). 

3.2. APPLICATIONS OF MEAN-FIELD THEORY WITH A SINGLE SOURCE 
OF HETEROGENEITY 

3.2. 1. Numerical bifurcation analysis using MFII 

As shown in Figure 3 the CA3 model network a makes transition 
from tonic firing to bursting as (i a pp) is varied. Similar transitions 
occur when g syn is varied. In this section, we use numerical bifur- 
cation analysis of MFII to determine the bifurcations involved in 
this transition, and the curves where they occur in the (4pp)-&yn 
parameter space. Since the mean-field system (Equations 50-54) 



consists of switching ODEs, this involves bifurcations of non- 
smooth systems as well as standard (smooth) bifurcations. A 
review of the theory of non-smooth systems can be found in 
di Bernardo et al. (2008). For the standard (smooth) bifurcations, 
the numerical bifurcation analysis is done in MATLAB (MATLAB, 
2012) using the MATCONT package (Dhooge et al, 2003). While 
it is possible to apply typical numerical continuation techniques 
to non-smooth bifurcations, primarily by defining alternate sets 
of test functions, and functions defining non-smooth limit cycles 
and equilibria for continuation (see Kuznetsov et al., 2003), 
implementation of these algorithms is outside of the scope of 
this paper. We opted instead to determine the non-smooth bifur- 
cations via direct numerical simulations. We compare the mean 
field theory results to those for the homogeneous system and to 
direct simulations of large networks. 

In Nicola and Campbell (2013) we carried out a numerical 
bifurcation analysis for a homogeneous network. The mean-field 
equations in this case, which are the same as MFI with (/app) 
replaced by I app , indicate that the transition from tonic firing to 
bursting occurs via the following sequence of bifurcations. The 
stable bursting limit cycle is created in a saddle node bifurcation 
of (non-smooth) limit cycles. The smaller, unstable limit cycle 
becomes smooth in a grazing bifurcation and then disappears in 
a subcritical-Hopf bifurcation which destabilizes the equilibrium 
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FIGURE 6 | Comparison between the bifurcation structure of 
homogeneous and heterogeneous networks using mean-field models. 

The parameters are as in Table 1, except the applied current which is 
normally distributed with mean and variance as follows and g S yn which varies 
as shown. (A) MFI, / app = 2000 pA. (B) MFII, (/ app ) = 2000 pA, a, = 500 pA. 
(C) MFI, low g syn bifurcation sequence (3D view). (D) MFI, low g syn 
bifurcation sequence, (2D view). In (A,B), the curved blue lines denote the 
value of the equilibrium point, which corresponds to tonic firing in the 
network. The vertical black/blue lines denote the amplitude range for the 
stable/unstable limit cycles, respectively. These correspond to bursting if the 
amplitude reaches zero, otherwise they correspond to the network having an 
oscillatory average firing rate. (A) Homogeneous case. Numerical bifurcation 
analysis of MFI displays two subcritical Hopf bifurcations: one at a low g sy n 
value and one at a high value. (B) Heterogeneous case. Numerical bifurcation 
analysis of MFII also displays two Hopf bifurcations, but the one at the low 



g syn value is supercritical. This makes bursting at low g syn values less robust 
in the heterogeneous case as discussed in the text. The non-smooth 
sequence of bifurcations for the homogeneous network is expanded upon in 
(CD). In (C), we continue the smooth unstable limit cycle (blue) from the 
Hopf bifurcation at low g syn values until the continuation halts (red). This 
occurs close to the grazing bifurcation with the switching manifold. We use 
direct numerical simulations to continue the stable non-smooth limit cycle 
(green). Key points in this bifurcation sequence are shown in the phase plane 
in (D). A smooth unstable limit cycle expands and grazes the switching 
manifold at approximately g syn = 52.71 nS, and persists to collide in a 
non-smooth saddle-node of limit cycles at approximately g syn = 55.11 nS. 
This sequence of bifurcations happens very rapidly in the parameter space, 
leaving a narrow region of bistability between the tonic firing and bursting 
solutions. Note that the totally unstable-non smooth limit cycles in (D) are 
computed via direct simulation of the time reversed MFII system. 



point corresponding to tonic firing. This transition is shown in 
Figure 6A when I is held fixed and g syn is varied. The bursting 
limit cycles are created at a low g syn value and destroyed at a high 
g syn value. Details of the sequence of bifurcations for lowg 5yn are 
illustrated in Figures 6C,D. The sequence for high g syn is the same 
but reversed. 

Using MFII, we numerically confirm that, as for the homo- 
geneous network, the mean-field system of the heterogeneous 
network undergoes a Hopf bifurcation as the network transitions 
from tonic firing to bursting. However, as shown in Figure 6B, 
with (/app) held fixed the transitions for low g syn and high g syn 
are not the same. For high g syn the transition is the same as the 
homogeneous case. For low g syn the transition occurs via the fol- 
lowing sequence of bifurcations. A supercritical Hopf bifurcation 
destabilizes the equilibrium point corresponding to tonic firing 
and creates a stable limit cycle. This limit cycle is smooth and 



hence corresponds not to bursting, but to firing with an oscil- 
latory firing rate. This limit cycle then grows until it becomes 
a non-smooth, bursting limit cycle in a grazing bifurcation (in 
Figure 6B this occurs at g syn ~ 150). We verified this prediction 
of the mean-field model by running direct simulations of a net- 
work of 10,000 neurons with fixed (I a pp), crj while varying the 
gsyn value. As shown in Figure 7A, when the steady state mean 
variables are plotted vs g syn , the supercritical nature of the Hopf 
bifurcation in the large network is apparent. 

To further investigate the heterogeneous case, we used 
MATCONT to carry out two parameter continuation of the Hopf 
bifurcation for MFII with four different values for the stan- 
dard deviation of J app : oi = 250,500,750, and 1000 pA. As shown 
in Figure 7B, in all cases there appears to be a codimension-2 
Bautin (or generalized Hopf) bifurcation on the two-parameter 
Hopf bifurcation curve, with the Hopf being supercritical on the 
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FIGURE 7 | Comparison between numerical bifurcation analysis of MFII 
and direct simulation of the full network. The parameters are as in 
Table 1, except g syn varies as discussed below and the applied current which 
is normally distributed with mean and variance as discussed below. (A) 
Simulations of a network of 10,000 neurons with {/ ap p) and cr/ as shown were 
run at discrete values of g syn for 2000 ms. The last 400 (ms) of simulation 
time is plotted (in red), showing the stable limit cycle oscillation for different 
g S yn values. This is compared to numerical continuation of the MFII limit 
cycle and equilibrium (in green and blue). Both the actual network and MFII 
appear to undergo a supercritical Hopf bifurcation for low g values and a 



subcritical Hopf for high g values. (B) The Hopf bifurcation curves for the 
mean-field systems with a/ as shown. Red denotes supercritical Hopf 
bifurcations and blue denote subcritical Hopf bifurcations. The black circles 
denote codimension 2 Bautin bifurcation points. (C) Simulations of a network 
of 1000 neurons run on a discrete mesh of </ a pp) and g syn values. The 0% 
(dotted line) and 100% (solid line) network bursting contours for a ; = 0, 250, 
500, 750, and 1000 pA are colored in black, magenta, blue, green, and red, 
respectively. The curves are spline fits to the actual contours. (D) MFII Hopf 
bifurcation curves and spline fits to the 0% bursting and 100% bursting 
contours of the actual network for cr/ = 1000. 



left boundary before this point and subcritical after. By con- 
trast, in the homogeneous case (07 = 0 line in Figure 7B) the 
bifurcation is subcritical everywhere on the two-parameter Hopf 
curve. 

Further verification of the mean-field results can be found in 
the direct numerical simulations of a network of 500 neurons 
shown in Figure 7C. The simulations were run on a 50 x 50 mesh 
in the g syn vs (J a pp) parameter space, using five different values for 
the standard deviation of I app : cr/ = 0, 250,500,750, and 1000 pA. 
Note that cr/ = 0 is the homogeneous network. The proportion 
of bursting neurons, pburst, was computed using Equation (A6) 
(see Appendix B). The 0 and 100% bursting contours can be seen 
in Figures 7C,D. In these figures, there appear to be two kinds 
of transitions from tonic firing to bursting. Along the (lower) left 
part of the boundary of the bursting region, the transition is grad- 
ual: the proportion of bursting neurons gradually increases from 
0 to 100%. Along the rest of the boundary, however, the whole 
network transitions to bursting simultaneously. This agrees with 
the prediction from the mean-field model that two different bifur- 
cations occur along the bursting boundary. Note also that the 



size of the entire bursting region and the 100% network bursting 
region get smaller as the level of heterogeneity (07) increases. 

Let us reiterate the primary differences between supercritical 
and subcritical Hopf induced bursting seen in Figures 6, 7. First, 
the subcritical case allows for bursting a lower g syn values than the 
supercritical case. This is because in the subcritical case bursting 
is initiated via a a saddle-node of limit cycles bifurcation which 
occurs to the left of the Hopf bifurcation, while in the supercritical 
case, bursting starts to the right of the Hopf in a grazing bifurca- 
tion. Second, the transition to bursting is sharp in the subcritical 
case and gradual in the supercritical case. The supercritical Hopf 
bifurcation is consistent with the gradual transition from burst- 
ing to firing seen in Figures 7C,D. When only a few neurons are 
bursting and the rest have oscillatory firing rates the correspond- 
ing mean behavior is a limit cycle with small amplitude. As more 
and more neurons become bursting this increases the amplitude 
of the limit cycle of the mean behavior until it grazes the switch- 
ing manifold. In the subcritical case, the saddle-node of limit 
cycles involves large amplitude limit cycles, corresponding to all 
the neurons being in the bursting state. 
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The bifurcation curves in Figure 7B are both qualitatively and 
quantitatively accurate descriptions of the behavior of the actual 
network. For example, in the actual network simulation, the 
bursting region decreases as 07 increases (see Figure 7C). This 
same behavior is displayed by the MFII equations, albeit to a 
greater degree, as shown in Figure 7B. However, there is a greater 
degree of quantitative error for lower values of g syn and larger val- 
ues of 0[. In particular, for a fixed value of cr; MFII predicts that 
the Hopf bifurcation occurs at a higher value of g syn than occurs 
in the real network (compare Figures 7B,D directly) and this pre- 
diction error seems to increase as 01 increases. This is why MFII 
indicates the network should be tonically firing when ai is high 
(in Figure 3D, for example). 

Taken together, these results indicate that for small g syn , net- 
work induced bursting via adaptation is not robust to hetero- 
geneity in the applied current. This occurs for qualitative and 
quantitative reasons, both related to the Hopf bifurcation associ- 
ated with the left boundary of the bursting region. Qualitatively, 
the addition of heterogeneity causes this bifurcation to change 
from subcritical to supercritical making the bursting less robust 
for small g syn values. Quantitatively, the g syn value of this bifurca- 
tion increases when the heterogeneity becomes stronger, while the 
value of the Hopf bifurcation associated with the right boundary 
does not change appreciably. Thus the size of the bursting region 
decreases with increasing heterogeneity. 

3.2.2. Bifurcation types and manifolds using MFII I 

It is difficult to use MATCONT with MFIII as MFIII is an infi- 
nite dimensional dynamical system, as it is a PDE. However, the 
existence of equilibrium points can be determined using standard 
root finding algorithms. While direct bifurcation analysis is diffi- 
cult to implement in this situation, one can use properties of the 
firing rate to describe, qualitatively and quantitatively, any tran- 
sitions between network states. This will be the approach of this 
section. To begin, we consider networks that are tonically firing, 
we then proceed to the study of bursting networks. 

For a network of neurons with heterogeneity in the parame- 
ters, even if all the neurons are tonically firing, one cannot find 
a steady state firing rate for the network, as in the case of a 
homogeneous network. The parameter heterogeneity creates a 
distribution of steady state firing rates across the network. While 
the mean-field equations by themselves can only determine the 
mean of this distribution, with an added assumption we can 
approximate the distribution of steady state firing rates for the 
network with a great degree of accuracy. 

Consider a network with just one heterogeneous parameter, p. 
Assume that the steady state firing rate of each neuron in the net- 
work can be related to its value for the heterogeneous parameter: 
R ; = g(|3). Assume further that one can approximate this function 
by the steady state value of (R,(t) |P> : 



Pr{t), through the standard theorem on transforming random 
variables: 



(60) 



This is easily determined through direct simulation of MFIII, 
Equations (56-58), until the system reaches steady state. Treating 
g as the transformation of a random variable, one can deter- 
mine the steady state distribution of firing rates in the network, 



Pn(r) = pp(g (r)) 



dr* 



V) 



(61) 



which can be found in any standard textbook on probability the- 
ory (such as Renyi, 1970). Note that we must assume that (-R,|P) 
is monotonic and invertible for this procedure to be valid. 

We carried out this computation for a network of 1000 neu- 
rons with a normal distribution in either I, g, or Wj ump . Details 
of the implementation can be found in Appendix D. We numer- 
ically determined the distribution of steady state firing rates for 
the neurons in the full network through 



1 

R, = , i = 1,2, 



. N 



(62) 



where ISIy as t is the last interspike interval for the ith neuron mea- 
sured from a lengthy (1000 ms) simulation. Figure 8 shows the 
results of the two approaches. The blue curve in the left column 
shows the distribution of parameter values. This is used in MFIII 
to calculate the predicted distribution of firing rates, which is 
the dashed red curve in the right column. The solid blue curve 
in the right column is the computed distribution of firing rates 
from numerical simulation of the full network equations. There 
is excellent agreement between the firing rate distributions in all 
cases. 

We carried out the same computations with a bi-modal distri- 
bution in I, g, or Wj ump , generated by mixing normal unimodal 
distributions (see Appendix C). This is one way of representing 
a network with two subpopulations of neurons with different 
parameters. The mean field approach again gives an excellent 
approximation to the qualitative and quantitative properties of 
the steady state distribution of firing rates, as shown in the right 
column of Figure 9. 

The above approach is only valid when the network is tonically 
firing. In this situation the steady firing rates of the network and 
the individual neurons are constant. When the network leaves the 
tonic firing regime, however, these steady state firing rates become 
oscillatory. In the case of bursting, the amplitude of the oscilla- 
tion is large enough that the firing rate goes to zero for intervals 
of time. Whether or not the neurons are bursting, oscillatory fir- 
ing rates cannot be represented as a simple distribution of firing 
rates. However, with some additional work we can use the tools 
developed above to determine what proportion of neurons in 
the network is bursting. This is a statistical alternative to direct 
bifurcation analysis. 

From simulations of the full network, we know that when 
the variance in the heterogeneity is large enough, not all of 
the neurons necessarily display the behavior predicted from the 
mean-field equations. For example, Figure 10 shows simulations 
of a network where the mean-field equations display an oscilla- 
tory firing rate which does not quite go to zero. The spike time 
raster plot of the full network (Figure 10A) shows that some 
neurons are bursting while others are tonically firing with an 
oscillatory firing rate. However, simulation of the corresponding 
mean-field equations (Figure 10B, dashed line) reveals only an 
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FIGURE 8 | Heterogeneity in / app , g syn , or Wj ump leads to heterogeneity 
in the firing rate. As described in section 3.2.2, given the parameter 
distribution (solid curves, left column) MFIII can be used to estimate the 
corresponding distribution of firing rates in the network (dashed curves, 
right column). As described in section 3.2.3, given the steady state firing 
rate distribution of an actual network (solid curves, right column) MFIII can 
be used to estimate the parameter distribution in the network (dashed 
curves, left column). The network firing rate distribution is estimated using 
a histogram. The calculations were carried out on a network of 1000 
neurons. Parameters, other than those given below, can be found in 
Table 1. (A) Distribution of / app with </ app ) = 4500 pA, ct, = 1000 pA. (B) 
Distribution of g S yn with (g S yn) = 200 nS, a g = 50 nS. (C) Distribution of 
l/V, ump with (l/Vjump) = 200 pA, o w = 50 nS. 
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FIGURE 9 | Bimodal distributions in / app , g syn , and M/j ump lead to 
bimodal distributions in the firing rate. These bimodal parameter 
distributions are generated through distribution mixing of two normal 
subpopulations with standard deviations and means as given below. See 
Appendix C for details. The distribution of the firing rate or the distribution 
of the parameter can be computed using MFIII if one knows the 
complementary distribution. See sections 3.2.2 and 3.2.3 for details. Curve 
descriptions are as given in Figure 8. Parameters, other than those given 
below, can be found in Table 1. The calculations were carried out on a 
network of 1000 neurons. (A) Distribution of / app with m = 4500 pA, 
at = 500 pA, H2 = 7000 pA, a 2 = 1000 pA, m = 0.7. (B) Distribution of g sy n 
with fit = 100 nS, at = 30 nS, n 2 = 300 nS, o 2 = 50 nS, m = 0.4. (C) 
Distribution of l/Vj ump with [x q = 300 pA, a-\ = 50 pA, H2 = 75 pA, 
CT 2 = 20 pA, m = 0.6. 



oscillation, not bursting. While this is consistent with the behav- 
ior of the network mean variables (Figure 10B, solid line), we 
have lost the information that some of the neurons are bursting. 
Similarly, one can find examples where the mean-field equations 



exhibit bursting, but not all neurons in the network are burst- 
ing. Thus, it would be useful to have more information about 
individual neuron behavior. MFIII can be used to obtain such 
information. 
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FIGURE 10 | Visualizing a limit cycle in a heterogeneous network. 

Numerical simulation of MFIII and a network of 1000 neurons with 
heterogeneity in the applied current. Parameters are as given in Table 1 
except g syn = 200 nS, (4 PP ) = 1000 pA and a, = 4400 pA. (A) Raster plot 
for 50 randomly selected neurons of the network arranged in order of 
increasing current. Some of the neurons are bursting, while others are 
tonically firing, albeit with an oscillatory firing rate. (B) In the mean 
variables, the steady state behavior of both the network and MFIII is an 



oscillation. (C) As MFIII is a partial differential equation, the steady state 
"limit cycle" is actually a manifold of limit cycles, foliated by the 
heterogeneous parameter p = / app . Part of the manifold has cycles with 
(R|P) = 0 for an extended period of time (in blue). The other part contains 
limit cycles that have (fl|fi) ^ 0 during the entire oscillation. We can 
classify neurons with the parameter values in blue as bursting, and those 
in green as oscillatory firing. (D) Averaging the limit cycle in (C) with 
respect to p yields the mean limit cycle. 



In the situation described above, the steady state mean network 
firing rate is oscillatory. We will denote this oscillatory solution as 
y. We interpret y to be the limit cycle parameterized by y(f) with 
(s(y(f)), (w(y(O)IP)) being the graph of the limit cycle in phase 
space. We will denote the period of the limit cycle as T. In this 
case, the "steady state" firing rate for neuron i will be a periodic 
function of time: R,(t), which depends on y and the value of the 
parameter p associated with the neuron: R,(t) = g(fi, y(f))> for 
t € [0, T]. To proceed, we make the same assumption as above, 
that 

g(p\y(f))~<£,(Y(f))IP> (63) 

where (^i(y(f))IP) is the oscillatory firing rate associated with 
the steady state limit cycle y in MFIII. An example of the graph 
of the steady state limit cycle derived from MFIII is shown in 
Figure 10C. In this visualization we can clearly see that part of the 
network is bursting (blue) while the rest is tonically firing with 
an oscillatory firing rate (green). Integration of this limit cycle 
over the heterogeneous parameter returns the "mean" limit cycle 
(Figure 10D). 

We now use this setup to approximate pbursb the proportion of 
neurons in the network that are bursting during the network level 
oscillation y. Noting that (-R,(y(0IP)) > 0, for all t 6 [0, T], the 



tonically firing neurons correspond to those P values for which 
CR;(Y(t))IP) > 0, for all t e [0, T\. Thus we define, p tonic , the 
proportion of tonically firing neurons in the network via 



Ptonic 



M 



min (R;(Y(t))IP> 
fe[0,T] 



0 p p (P)dp (64) 



where X is the usual indicator function. Similarly, the proportion 
of quiescent (non-firing) neurons is given by 



Pi 



H 



max (Hi(y(O)IP) 
te[0,T\ 



0 Pp (P)dp. 



(65) 



Recall that the bursting neurons correspond to those p values 
such that (_R(y(f))|P) = 0, for some subinterval of [0, X]. Thus 
we must have 

Pburst = 1 — Pq — Ptonic- (66) 

We compute these values as follows. First we numerically inte- 
grated MFIII until the steady state oscillation y is reached. This is 
an oscillation of (w(f)IP) and s(t). The corresponding oscillatory 
firing rate (P,(y(f))IP) is computed through Equation (58) as a 
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function of p on the limit cycle y. One can then determine 



m(p) = min <R(y(f))IP>, 
M(P)= max <R(y(f))|p>, 

K[0,T] 



and the integrals simplify to 



Ptonic= f h(«(P))/S|j(P)<* 

p q = [ fe(-M(P))pp(p)dp 



(67) 
(68) 

(69) 
(70) 



where h is the Heaviside function, with h(0) = 1. 

Numerical results for a network with single heterogeneous 
parameter are shown in Figure 11. For Figure 11 A, the pro- 
portion of bursting neurons was computed, using the method 
described above, at each point in a mesh on the g syn vs. (4pp) 
parameter space. This data was then used to generate the pburst 
contours. For Figure 1 IB the actual network was simulated to 
steady state at each point of a slightly coarser mesh. The propor- 
tion of bursting neurons at each point was computed according 
to Equation (A6) in Appendix B and used to generate the pburst 
contours. The results of the mean-field computation are both 
qualitatively and quantitatively accurate. In particular, MFIII 
recovers the gradual transition to bursting on the left boundary 
of the bursting region and the abrupt transition to bursting on 
the right. It should be noted that it is much faster, by approxi- 
mately an order of magnitude, to run a mesh of integrations over 
MFIII then it is to run mesh over an actual network. 

3.2.3. Inverting a steady state firing distribution to determine the 
distribution of parameters using MFIII 

Many parameters for neuron models are difficult to measure 
directly using electrophysiology. However, a distribution of fir- 
ing rates across a network of neurons is relatively easy to measure 
using intracellular recordings, or can be estimated using mea- 
surements from multi-electrode recordings and spike sorting 
algorithms, among other methods (Buzsaki, 2004; Grewe et al., 
2010). We have seen in the previous section that, given a dis- 
tribution of heterogeneities, MFIII can predict the steady state 
distribution of firing rates. Here we show that one can invert this 
process to yield a distribution of parameters given a steady state 
distribution of firing rates. 

We assume that only the firing rate distribution is known, and 
denote it Pr(t) as above. We then proceed as in the previous 
section, assuming that the steady state firing rate for a partic- 
ular neuron is some function of the heterogeneous parameters 
Rj = g(P) and that this function is well approximated by (RjIP)- 
Under these assumptions, one can solve for the distribution of 
parameters P using 



/op(P) = Pfi(g(P)) ^gS(P) 



(71) 



need to assume that (-R,|P) is differentiable for this procedure to 
be valid. 

The primary problem we face in using this approach to 
approximate the distribution pp(P) is that we need to determine 
the steady state values of the function (i?;|p). However, a cursory 
look at the equations for MFIII shows that these in fact depend 
on pp(P), the function we are trying to find, through the equation 
for s: 



~r~ Sjump 



f <Sj(t)IP)flp(P)dp. 



(72) 



Fortunately, however, this problem disappears when we look at 
the steady state value for s: 
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which follows from standard statistical theorems on the 
transformations of random variables (Renyi, 1970). Note that we 
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FIGURE 11 | The proportion of bursting neurons, Pburst f° r MFIII and an 
actual network. (A) Using the techniques outlined in the text, MFIII is 
used to compute the proportion of bursting neurons, Pb ur st at each point in 
a mesh over the parameter space. (B) Numerical simulations of a network 
of 500 neurons are used to compute the proportion of bursting neurons, 
Pburst- All the parameters for both the MFIII system and the actual network 
are identical (see Table 1), except that MFIII is run over a finer mesh. The 
results using MFIII are both qualitatively and quantitatively accurate. (A) 
MFIII, 07 = 500 pA. (B) Network, 07 = 500 pA. 
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f (RilP)/Op(P)dP = t s s jump (J?,). (73) 

h 



Here (R{) is the unconditioned steady state mean of the fir- 
ing rate distribution. This information is readily available, as we 
have assumed we know the steady state distribution, Pr(t), and 
determining the first moment is numerically trivial. 

Putting the expression for s into the steady state equation for 
(w|P) yields a set of coupled equations: 



(W|P) = TwWjump^ilP), 

[A 

0 



(74) 

< (wiw.p,]"' :H «HP),5,P)>o 

(75) 



H((w|P),s, P) <0 



These maybe solved for (w|P) and (i?,-|P) by discretizing in p and 
numerically solving the resulting system at each grid point with 
any standard root finding algorithm. 

Alternatively, one can set s to its equilibrium value in MFIII 
and numerically integrate the resulting equation: 



(w|P)' = -fl(w|P)+w jump (i?i«IP) 



(76) 



<«i(0 IP) 



0 



> ( : lll .. T , %lllt , W) .MP), P ) ]" :H(HP),"5,P)>0 

(77) 

H((w|P),s,P) <0 



until it reaches steady state, which will determine (w|P) and 
(.RilP). Note that this approach will only work if the tonic 
firing equilibrium of the original mean-field system MFIII is 
asymptotically stable. 

We have implemented this approach as follows. A network of 
1000 neurons is numerically integrated until it reaches its steady 
state firing rate. The distribution of firing rates over the network is 
found as described in the previous section. The density function 
for this distribution, pj?(r), is then estimated using the firing rate 
histogram. Equations (76), (77) are numerically integrated until 
they reach steady state. We then substitute the estimate of pr (r) 
and the approximation (R,|P) ofg(P) into (71) to determine the 
parameter distribution pp(P). See Appendix D for more details. 
Our results for unimodal and bimodal distributions are shown 
in Figures 8, 9, respectively. In the right column of each figure, 
the solid blue curve is the distribution of steady state firing rates 
from integration of the full network. In the left column of each 
figure the dashed red curve is the estimate of pp(P) found using 
the procedure above, while the blue curve is the actual param- 
eter distribution used in the network simulation. We note that 
no information about the distribution of parameters is known 
in the estimation procedure, yet the numerical results are very 
accurate in both the unimodal (Figure 8) and the bi-modal case 
(Figure 9). 

Perhaps most interesting is that we can extend this tech- 
nique to estimate the individual neuron parameter values, p,, i = 
1, . . . , N. This again follows from the assumption that g(P) = 
(Ri | P) is the function that transforms the random variables P; into 
Ri. If the function is invertible, then we can compute the individ- 
ual p, through numerically inverting the steady state (R;|P). For 



example, when this technique is applied to a network where the 
only source of heterogeneity is J, the mean relative absolute error 
in the predicted values I, versus the actual values J; is only 0.6%. 
The details about how to numerically invert for the individual 
parameter values are included in Appendix D. 

While network level inversion of a single heterogeneous 
parameter is an important step forward, this is performed under 
very strong assumptions. In particular, when performing this 
inversion, all of the heterogeneity in the firing rates is assumed 
to come from a single parameter. Additionally, all the other 
parameters are assumed to be known. These two assumptions are 
exceptionally strong and one has to take great care in inverting 
actual recorded firing rates from neurons that they be reasonably 
satisfied. 

3.3. MEAN-FIELD APPLICATIONS WITH MULTIPLE SOURCES OF 
HETEROGENEITY 

In order for the mean-field applications to be useful for realis- 
tic neuronal networks, one needs to consider heterogeneity in 
more than one parameter. Recall that the mean field systems 
derived in section 2.3 are valid for multiple heterogeneous param- 
eters, one simply considers p to be a vector instead of scalar. This 
presents some difficulties in implementation which we discuss in 
the section. The examples we consider will have 2 or 3 sources of 
heterogeneity, primarily in the parameters J, g and d. 

Recall that MFII is given by the Equations (50-54). The main 
difficulty in dealing with MFII lies with the integral terms, which 
are now multiple integrals. For example: 

(Ri)= / / ...[ <a i |p)pp(p 1 ,p 2 ,...B p )dp p ...dp 2 p 1 (78) 

where p is the number of heterogeneous parameters. In order to 
numerically integrate or carry out bifurcation analysis on MFII, 
these multiple integrals must be evaluated. We have found that 
this is most easily done using a Monte-Carlo numerical inte- 
gration scheme. Once this is implemented, bifurcation diagrams 
can be generated exactly as for the case of one parameter het- 
erogeneity: the equilibrium points and smooth limit cycles are 
continued using MATCONT, while the non-smooth limit cycles 
are generated using numerical simulations. 

The integral term in MFIII can be dealt with in a similar way as 
to that for MFII. Once this is implemented, the steady states and 
network properties can be determined as described in section 3.2, 
while numerical simulations can be used to follow stable periodic 
solutions. We will use this approach later on in our case study on 
adaptation induced bursting. 

Mean-field III can also be used to determine steady state fir- 
ing rates following the procedure in section 3.2.2, however, one 
now has to discretize the equations over a multi-dimensional 
mesh. While this approach is feasible, we found it is more efficient 
to predict the steady state firing rates of the individual neurons 
through the following interpolation scheme. Given knowledge of 
the parameter distribution, we generate sample points ft ; from 
this distribution. We then generate a steady state firing rate for 
each sample point, _R, = g(P,), where g is determined from MFIII 
as described in section 3.2.2. We interpolate over the (p\, i?,) 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



December 2013 | Volume 7 | Article 184 | 19 



Nicola and Campbell 



Mean-field models for heterogeneous networks 



ordered pairs to determine the firing rates of the individual neu- 
rons, given knowledge of their parameter values. If we only need 
the distribution of the firing rates, then the distribution of the 
_R, is an estimate of this, without need for interpolation. This 
approach has been applied to a network of 1000 Izhikevich neu- 
rons with three simultaneous sources of heterogeneity, as shown 
in Figure 12. 

The one application we found difficult to extend to the case 
of multiple sources of heterogeneity was the mapping of the 
distribution of steady state firing rates to the distribution of 
parameters. There is a fundamental difficulty with this inversion 
problem: the firing rate distribution is one-dimensional, but the 



distribution of parameters is multidimensional. Thus, we leave 
further investigation of this problem for future work. 

3.3.1. Bifurcation analysis with multiple sources of 
heterogeneity — case study 

To conclude our work, we consider a realistic model for a CA3 
hippocampal network of pyramidal cells. Hemond et al. (2008) 
classify CA3 pyramidal cells into three types: weakly adapting, 
strongly adapting and intrinsically bursting. We will focus on the 
effect on network bursting of having two subpopulations: one 
strongly adapting and one weakly adapting. We use the Izhikevich 
model (Equations 4-6) with the parameters set up by Dur-e- 
Ahmad et al. (2012) (see Table 1), but include heterogeneity in 
^app> gsyn and the adaptation parameters Wj ump , xw- The param- 
eter distributions are generated through distribution mixing (see 
Appendix B) of normal distributions with the parameters given 
in Table 2. We have treated the mean values of 7 a pp and g syn from 
the strongly adapting subpopulation as the bifurcation parame- 
ters. We also varied the proportion of strongly adapting neurons 
in the population, i.e., parameter p in Equation (A8). 

The 0 and 100% bursting contours for simulations over the 
two parameter mesh in the (-f a pp), (gsyn) are shown in Figure 13 
for both the full network (Figure 13A) and MFIII (Figure 13B). 
Numerical bifurcation analysis of MFII (not shown) confirms 
that the bifurcations are similar to when J app was the only source 
of heterogeneity, in particular, on the left boundary of the burst- 
ing region the Hopf bifurcation is supercritical, while on the right 
it is subcritical. 

As shown in Figure 13, when the proportion of strongly adapt- 
ing neurons is decreased, the bursting region decreases. However, 
unlike previous results (Nicola and Campbell, 2013, Figure 10), 
the decrease seems to be more pronounced in the high g syn region. 
This is likely due to having truly heterogeneous distributions 
of parameters, as opposed to splitting a network into two dif- 
ferent homogeneous subpopulations as was done in Nicola and 
Campbell (2013). In all cases, it appears that heterogeneity shifts 
the bursting region to higher values of g syn , outside the range 
of biologically plausible conductances described in our previous 
work (Nicola and Campbell, 2013). 

4. DISCUSSION 

Building on the mean-field framework for networks of homo- 
geneous oscillators, we extended the mean-field approach to 



Table 2 | Table of parameters for the strongly and weakly adapting 
heterogeneous subpopulations. 
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FIGURE 12 | A bimodal distribution in / app together with unimodal 
distributions in <y S yn and Wj U mp as shown in (A) yields a unimodal 
distribution in the firing rate [dashed curve in (B)]. The mean-field 
equations give a good estimate of this distribution [solid curve in (B)]. 
Simulations are for network of 1000 neurons, parameter given in Table 1 
except /, g, and l/Vj U m P - These are distributed with o w = 50, {W lump ) = 200, 
o g = 50, <g syn ) = 200, m = 0.5, (/ app ,i) = 7500, a/,-, = 200, </ apPi2 ) = 5500, 
a/ 2 = 500. Other combinations of bimodal and unimodal parameter 
distributions may yield bimodal firing rate distributions. Note that this is 
different than the situation shown in Figure 9, where a single bimodal 
source of heterogeneity yielded a bimodal firing rate distribution. 
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FIGURE 13 | Case study: the proportion of bursting neurons in a 
network with both strongly adapting and weakly adapting neurons. 

Parameters are as in Table 1 except the / app , g syn , IA^ ump , and tw which 
have bimodal distributions generated by distribution mixing with 
parameters given in Table 2. The parameter psa represents the proportion 
of strongly adapting neurons in the network (see Equation A8). The dashed 
line is the 0% bursting contour, while the dotted line is the 100% bursting 
contour. As shown in both the network (A) and MFIII (B), the bursting 
regions becomes significantly smaller when the proportion of strongly 
adapting neurons decreases. In all cases, the curves are smoother spline 
fits to the actual contours. 



networks of heterogeneous oscillators. This was accomplished 
through the derivation of three separate mean-field systems, 
MFI, MFII, and MFIII, with differing applications and regions 
of validity. We successfully applied numerical bifurcation anal- 
ysis to MFI and MFII to aid in the understanding of the dif- 
ferent behaviors that heterogeneous networks can display, and 
how they transition between these different types of behav- 
iors. More importantly, however, we have surpassed the natural 
limitation of mean-field systems: that they can only provide 
information about the first moments. With a few additional 
tools, we used MFIII to derive information about distributions 



of firing rates, and even parameters, given some basic 
knowledge. 

Other researchers (Hansel and Mato, 2003; Vladimirski et al., 
2008; Hermann and Touboul, 2012) have derived firing rate dis- 
tributions for heterogeneous networks, however, these have been 
derived under differing assumptions. For example, the heteroge- 
neous mean field systems studied by Hansel and Mato (Equations 
(5.5-5.7) in Hansel and Mato (2003)) and Hermann and Touboul 
(Equations (1, 2) in Hermann and Touboul (2012)) have similar 
integral terms to our MFII, however, they are firing rate models. 
Our models are current/conductance based models. The differ- 
ence between these two types of equations arises from which time 
scale is the fastest, that of the synaptic current, or the firing rate. 
If the firing rate time scale is assumed to be the fastest, then a dif- 
ferential equation for the synaptic current can be obtained (as in 
our case). If the time scale of the synaptic current is assumed to 
be the fastest, then one obtains firing rate equations, as in Hansel 
and Mato (2003) and Hermann and Touboul (2012). The fact 
that these two different limits result in different kinds of equa- 
tions were first highlighted in Dayan and Abbott (2001) (section 
7.2). Additionally, no adaptation is contained in the rate models 
in Hansel and Mato (2003) and Hermann and Touboul (2012). 
Finally, it is likely that the firing rate models shown in Hansel 
and Mato (2003) cannot display period doubling bifurcations as 
they have a similar structure to MFII, which misses out on the 
more complicated bifurcations of the actual network that MFIII 
can reproduce, due to its PDE nature. The firing rate models 
in Hermann and Touboul (2012) have a complicated bifurca- 
tion structure as they involve two subpopulations (excitatory and 
inhibitory) leading to a four dimensional ODE system. 

The model of Vladimirski et al. (2008) is formulated in terms 
of an input-output relation for the synaptic conductance, so has 
a different structure than ours. It involves a distribution of the 
synaptic depression variable so has some aspects similar to our 
MFIII, however, no PDE governing the evolution of this variable 
is derived. 

Dur-e-Ahmad et al. (2012) studied adaptation induced burst- 
ing in a network of homogeneous Izhikevich neurons, with 
parameters determined from experimental data on CA3 pyrami- 
dal neurons. They showed that, if the adaptation is strong enough, 
network bursting occurs in large regions of the parameter space 
consisting of the synaptic conductance, g syn , and the applied cur- 
rent, 7 app . In Nicola and Campbell (2013) we showed that the 
transition from tonic firing to bursting involves a saddle-node 
bifurcation of non-smooth limit cycles, followed by a grazing 
bifurcation and a subcritical Hopf bifurcation. For fixed 7 app 
greater than rheobase but sufficiently small, there is one transi- 
tion from tonic firing to bursting at a low g syn value and another 
from bursting back to tonic firing at a higher g syn value. Thus 
the bursting region is a closed semi-circular region in the g syn , 
7 app parameter space. In Nicola and Campbell (2013) we showed 
that the size of this bursting region is reduced if the network is 
split into two homogeneous subnetworks, one strongly adapting 
and one weakly adapting. Here, we used the tools we developed 
to investigate how this adaptation induced network bursting is 
affected by heterogeneity in the parameters. Somewhat surpris- 
ingly, we have found that adaptation induced network bursting 
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is not very robust to heterogeneity. This has been confirmed by 
direct simulations of the full network, bifurcation analysis using 
MFII and analysis of the proportion of bursting neurons in the 
network using MFIII. This lack of robustness is caused by two 
changes to the homogeneous case: 

1. The low g syn Hopf bifurcation point moves toward higher 
values, thereby decreasing the size of the bursting region. 

2. The low gsyn Hopf bifurcation switches from subcritical to 
supercritical. This has two effects: 

• The bifurcation direction changes, eliminating the bursting 
at conductance values less than the bifurcation value. 

• The initial limit cycles created by the bifurcation are small 
amplitude oscillations in the firing rate as opposed to full 
bursts, thus the transition to bursting moves to even higher 
conductance values. 

Further, in networks with both weakly and strongly adapting neu- 
rons, heterogeneity caused the high g syn Hopf bifurcation value 
to decrease when the proportion of strongly adapting neurons is 
reduced. 

Let us now put our results in the context of experimental 
results on the CA3 region. Bursting is often seen in these stud- 
ies (Andersen et al., 2006, section 5.3.5). When the neurons have 
their synaptic inputs blocked, however, the majority (~80%) of 
these pyramidal neurons do not display bursting, but different 
degrees of spike frequency adaptation (Hemond et al, 2008). 
Thus, it would seem that adaptation induced network bursting 
should play a role in the CA3 network. However, the biophysi- 
cally important part of the parameter region is in the low g syn 
region (Nicola and Campbell, 2013). When this fact is taken in 
conjunction with our results described above, this would seem to 
weaken the case that adaptation induced network bursting is the 
only source of bursting in CA3 networks. Some other mechanism 
seems necessary. 

In their study of hippocampal CA3 pyramidal neurons, 
Hemond and colleagues note that roughly 20% of pyramidal 
neurons were intrinsically bursting. That is, the neurons burst 
without any synaptic input for some input current. It may by pos- 
sible that a small subpopulation of intrinsically bursting neurons 
can facilitate bursting in the rest of the network, however, this 
would depend on the conductance values connecting this partic- 
ular subpopulation to the rest of the network. This hypothesis can 
be tested relatively easily using a mean-field approach. All that is 
required is to fit a two-dimensional adapting model to the intrin- 
sically bursting neurons. This is feasible, as has been previously 
noted, all the two-dimensional adapting neurons can be turned 
into intrinsically bursting neurons by simple parameter changes 
(Izhikevich, 2003). The conductance parameter connecting this 
subpopulation to the rest is best treated as a bifurcation param- 
eter, with some estimate of the range in which it lies in from 
physiological data. 

While an intrinsically bursting subpopulation is the most 
promising avenue of study with regards to hippocampal burst- 
ing, synaptic depression has also been shown to induce bursting 
in oscillators that cannot otherwise display this behavior. In a 



model of the developing chick spinal cord, Vladimirski et al. 
(2008) found that heterogeneity actually makes the bursting more 
robust, as opposed to less as we have found. Thus it is possible the 
synaptic depression induced bursting is more robust to hetero- 
geneity than adaptation induced bursting. However, in this study 
the heterogeneity was via a uniform distribution in the applied 
current (as opposed to the Gaussian distributions we consider) 
and typically (4 PP ) was close to rheobase, which could also be 
factors in their results. 

In addition to area CA3 in the hippocampus, adaptation 
induced bursting has also been suggested as a possible mechanism 
for the generation of velocity controlled oscillators ( VCO's) in the 
entorhinal cortex by Zilli and Hasselmo (2010). The VCO's burst 
at frequencies that vary with the velocity of the animal. When 
a subset of VCO's signals are linearly added to a readout neu- 
ron, an interference pattern emerges and a grid cell is formed. 
Zilli and Hasselmo (2010) use a recursively coupled network of 
homogeneous Izhikevich neurons with adaptation variables given 
by Wj ump = 100 and l/xw = 0.03, parameters values such that 
adaptation induced bursting can occur. The network acts a sin- 
gle velocity controlled oscillator with the burst frequency varying 
with the velocity of an animal. This is done by fixing the g syn 
parameter at a specific value and inverting the F(I) curve, where 
F is the frequency of bursts and I is the homogeneous applied 
current to each neuron. Grid cells can be generated by using 
multiple networks and linearly adding their output currents to a 
read-out neuron. This was done under uncorrelated noisy inputs 
arriving to each neuron. However, Zilli and Hasselmo (2010) 
state that synchrony in the noise (which can come from the 
animals velocity signal for example) coming to each VCO net- 
work can disrupt grid-cell formation. Here, we have shown that 
a heterogeneous network of oscillators can still maintain a net- 
work level oscillation rate, even if the individual neurons have 
different behaviors. The network level behaviors are predicted 
from the mean-field systems. Given the fact that the individual 
neurons are heterogeneous, any synchronized noise input into 
the individual neurons should become increasingly desynchro- 
nized by the differing responses of the individual neurons. As a 
network level oscillation exists and the heterogeneity will likely 
desynchronize any noise coming to the individual oscillators, this 
is a plausible means of generating velocity controlled oscilla- 
tors. We leave this particular application of mean-field theory for 
future work. 

In either of these applications, a mean-field system may 
yield valuable insights as to the mechanisms of bursting and 
the parameter regions where they occur. By carefully choos- 
ing the appropriate bifurcation parameters and accounting for 
the level of heterogeneity in the neurons in the network, one 
can determine the bifurcation types and behaviors neurons 
in these different networks display, in addition to estimates 
of the different distributions to yield insights about the real 
cells. 
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APPENDICES 

APPENDIX A: COMPUTING ( V) AND ( V\fi) 

When the quasi-steady state approximation is applicable, a very 
convenient method emerges for computing the moments of v. In 
particular, the quasi-steady state approximation not only yields 
the steady state flux, but it also yields the steady state density: 



P(y) 



K(w),s) 



F(y) — (w) + gs(e r — v) + 7 



where J((w) , s) is the steady state flux ((Rj(t)) for H((w) , s) > 0). 
Obviously with the density, we can compute the moments of any 
arbitrary function of v. For example, (v) is given by 



fv m 



J((w), s)vdv 



) - <w) + gs(e,- -v)+I 



if H((w),s) >0(A1) 



However, once we cross the switching manifold, Equation (Al) is 
no longer valid. An approximation has to be made to compute (v) 
in this region. In particular, we will assume that (F(v)) = F((v)) 
and that the dynamics of (v) are fast relative to s and (w) . If so, 
then the dynamics of (v) are approximately given by 

(v)' « F((v)) - (w) + gs(e r — (v)) +1 (A2) 

and we can solve the steady state equation for (v) 

F((v)) - (w) + gs(e r - (v)) + I = 0 

However, we have to be careful here. Based on the assumptions 
on F(v) we know that there are only two solutions to this system, 
and we need to solve for the stable one. For the Izhikevich neuron 
for example, this is given by 



(v_) 



a + gs 



J-H((w),s) 



(A3) 



Computation of the integral for the Izhikevich neuron, in con- 
junction with Equation (A3) yields the following equation for (v) 
over the entire (w) , s plane: 



(v) 



WO) 

2 



log 



(Vpeak-°^ £ ) 2 +H((w) 



■ V-H((w>, s) 



i) 



H((w),s) > 0 



: H({w),s) < 0 
(A4) 

The same formulas can actually be used for (v|P), with the inter- 
pretation that they are now explicit functions of P In particular, 
one can show that under tonic firing, (v|P) is given by the same 
formula as (v) when H((w), s, P) > 0. Additionally, the formula 
for (v) for H((w) , s, P) < 0 is also a satisfactory approximation 
for (v|P). To once again reiterate that while the formula for the 
mean (v|P) under heterogeneity is the same as the formula for (v) 
under homogeneity, this does not imply that (v) = (v|P). In fact, 
to compute (v), one has to compute the traditional integral for it: 



(v) 



: [ (v|P> /0 p(P)dp 



APPENDIX B: COMPUTING THE PROPORTION OF BURSTING NEURONS 
IN DIRECT SIMULATIONS 

Performing direct simulations of a network is fairly straight 
forward, but it is somewhat difficult to use the results of the sim- 
ulations to automatically classify a given neuron, let alone the 
entire network, as bursting or tonically firing. We considered var- 
ious classifiers and opted to use the ratio of the largest to the 
smallest interspike interval for a neuron when the network has 
reached steady state. For a bursting neuron, the ratio of its largest 
interspike interval (which is the interburst interval) to its small- 
est interspike interval should be large. Thus for the zth neuron we 
define 

maxISI 

A, = — — (A5) 
mm ISI 

where the max/min is determined after a suitable time period of 
steady state behavior (either bursting or tonically firing). We set 
a critical ratio, k c , and classify neuron i as bursting if X; is higher 
than this ratio. The critical ratio is typically taken to be 2. This is 
the single neuron classifier. We use the single neuron classifier to 
make an estimate for the total proportion of bursting neurons in 
the network as follows 



Pburst 



(A6) 



i= 1 



where h is the Heaviside function. We can use pburst to classify 
a network by picking some specific critical value of pburst as a 
threshold for a classifier, however, it is more informative to simply 
plot the contours of pburst for a set of simulations run over a mesh 
in the bifurcation parameters of interest. 

APPENDIX C: MULTIPLE SUBPOPULATIONS VIA MIXED 
DISTRIBUTIONS OF HETEROGENEITY 

In Nicola and Campbell (2013) we considered multiple sub- 
populations for heterogeneous networks by simulating discrete 
homogeneous subpopulations within a network. For example, we 
simulated two subpopulations with two different sets of param- 
eters corresponding to weakly adapting and strongly adapting 
neurons as a first attempt at studying inhomogeneous networks. 
However, a more realistic way of analyzing subpopulations in het- 
erogeneous networks, is through the use of mixed distributions. 
A mixed random variable Z is a function of two or more ran- 
dom variables. For example, if X and Y are random variables with 
probability density functions fx and/y, respectively, then 



X with probability p 
Y with probability 1 - 



(A7) 



is a mixed random variable with probability density function 

f Z (z)=pfx(z) + (l-p)f Y (z). (A8) 

Now consider a network where the adaptation jump size Wj um p 
comes from two different subpopulations. In one subpopulation, 
(wj um p) is large, and in the other, (wj U mp) is small. We can sim- 
ulate these two subpopulations within an individual (all-to-all) 
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Mean-field models for heterogeneous networks 



coupled network using heterogeneity and the density function 

/wjump(w) = pSAfsA(w) + (1 - p S A)fwA(w). (A9) 

Here p$A denotes the proportion of strongly adapting neurons, 
SA denotes the strongly adapting subpopulation and WA denotes 
the weakly adapting subpopulation. Note that these densities have 
higher moments than simply (wj ump > . It is the magnitude of these 
moments that determine whether or not the density in the het- 
erogeneous parameter we are considering is bimodal, indicative 
of multiple subpopulations. Using this approach we can analyze 
how the bimodality affects the steady state distribution properties. 
We can also use more of the data gained from real networks for 
the purposes of simulation, as opposed to arbitrarily classifying 
neurons as strongly adapting or weakly adapting. The parameters 
of the individual density functions can be approximated, along 
with p, using standard statistical approaches (maximum likeli- 
hood estimation, for example). More importantly, however, the 
same mean-field equations apply to a unimodal or a multimodal 
distribution of heterogeneity. 

APPENDIX D: DIFFERENTIATION AND NUMERICAL INVERSION OF 
PARAMETER DISTRIBUTIONS 

The procedure in section 3.2.2 requires not only the estimation of 
(-RilP), but also the calculation of its inverse and the derivative of 
this inverse. To calculate (-R,|P) the mean field equations (56-58) 
are discretized in p and the resulting equations are numerically 
integrated to steady state, which yields the steady state value of 
(R;IP;) at each mesh point, p ; -. The inverse g~ l (r) is calculated via 
numerically inverting the steady state (R,-|fy) as a function of Pj 
as follows. The MATLAB function interpl is used to interpolate 
values of P given values of R using the steady state ((R;|Pj), Py) 



mesh points. The derivative of the inverse is calculated using a 
finite-difference approximation over the mesh: 

dr \ r = r ~ rj+i-Tj 

This is then used to find pj?(r) at each mesh point via Equation 
(61). 

To implement the computations in section 3.2.3, Equations 
(76, 77) are discretized in p and numerically integrating to com- 
pute the steady state value of (R,-|P) at each mesh point, (_R, |P ; ). 
This is then used to compute ^(.R|B) through a first order 
finite-difference approximation over the discrete mesh: 

These two quantities are then used as approximations of g(P) 
and g'(P) in Equation (71) to find pp(P) at each mesh point. 

To estimate the parameter values for individual neurons, we 
take the discretized steady state firing rate, (R,\$j), calculated 
as indicated above. We then invert the functional relationship 
between (-R,|p\) and Py, and interpolate the p values of the indi- 
vidual neurons using their firing rate, i.e., we treat (-R;IPj) and Py 
as the (xj, yj) points to be interpolated. This yields an estimate 
of the parameter values for each individual neuron, unlike in the 
approach for section 3.2.2 which yielded an estimate of the overall 
distribution. Note that all that is required for this computation is 
knowledge of the steady state firing rate for each neuron. We use 
the griddata and the interpl functions in MATLAB to perform 
the interpolation. 
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